#!/usr/bin/env python
# -*- coding: utf-8 -*-

""" By Martin Senande-Rivera
    For Towards and atmosphere more favourable to firestorm development in Europe """

import os
import glob, sys
import numpy as np
import xarray as xr
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from mpl_toolkits.basemap import Basemap, maskoceans
import matplotlib.colors as mcolors
import datetime as dt
from matplotlib.colors import LinearSegmentedColormap
import cmaps
from pandas import Series, DataFrame
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from matplotlib.ticker import FormatStrFormatter


font = {'size'   : 5}
matplotlib.rc('font', **font)

########## FWI ##########
path_FWI_EC_EARTH='../4-Projections_firestorm-risk/1-FWI/COSMO-crCLIM-v1-1/EC-EARTH/'
path_FWI_HadGEM2_ES='../4-Projections_firestorm-risk/1-FWI/COSMO-crCLIM-v1-1/HadGEM2-ES/'
path_FWI_MPI_ESM_LR='../4-Projections_firestorm-risk/1-FWI/COSMO-crCLIM-v1-1/MPI-ESM-LR/'
path_FWI_NorESM1_M='../4-Projections_firestorm-risk/1-FWI/COSMO-crCLIM-v1-1/NorESM1-M/'
path_FWI_EC_EARTH_2='../4-Projections_firestorm-risk/1-FWI/RCA4/EC-EARTH/'
path_FWI_HadGEM2_ES_2='../4-Projections_firestorm-risk/1-FWI/RCA4/HadGEM2-ES/'
path_FWI_MPI_ESM_LR_2='../4-Projections_firestorm-risk/1-FWI/RCA4/MPI-ESM-LR/'
path_FWI_NorESM1_M_2='../4-Projections_firestorm-risk/1-FWI/RCA4/NorESM1-M/'

# Read files
FWI_01 = xr.open_dataset(path_FWI_EC_EARTH+'FWI_ndays_trend_rcp85.nc')
FWI_02 = xr.open_dataset(path_FWI_HadGEM2_ES+'FWI_ndays_trend_rcp85.nc')
FWI_03 = xr.open_dataset(path_FWI_MPI_ESM_LR+'FWI_ndays_trend_rcp85.nc')
FWI_04 = xr.open_dataset(path_FWI_NorESM1_M+'FWI_ndays_trend_rcp85.nc')
FWI_05 = xr.open_dataset(path_FWI_EC_EARTH_2+'FWI_ndays_trend_rcp85.nc')
FWI_06 = xr.open_dataset(path_FWI_HadGEM2_ES_2+'FWI_ndays_trend_rcp85.nc')
FWI_07 = xr.open_dataset(path_FWI_MPI_ESM_LR_2+'FWI_ndays_trend_rcp85.nc')
FWI_08 = xr.open_dataset(path_FWI_NorESM1_M_2+'FWI_ndays_trend_rcp85.nc')
FWIpvalues_01 = xr.open_dataset(path_FWI_EC_EARTH+'FWI_ndays_pvalues_rcp85.nc')
FWIpvalues_02 = xr.open_dataset(path_FWI_HadGEM2_ES+'FWI_ndays_pvalues_rcp85.nc')
FWIpvalues_03 = xr.open_dataset(path_FWI_MPI_ESM_LR+'FWI_ndays_pvalues_rcp85.nc')
FWIpvalues_04 = xr.open_dataset(path_FWI_NorESM1_M+'FWI_ndays_pvalues_rcp85.nc')
FWIpvalues_05 = xr.open_dataset(path_FWI_EC_EARTH_2+'FWI_ndays_pvalues_rcp85.nc')
FWIpvalues_06 = xr.open_dataset(path_FWI_HadGEM2_ES_2+'FWI_ndays_pvalues_rcp85.nc')
FWIpvalues_07 = xr.open_dataset(path_FWI_MPI_ESM_LR_2+'FWI_ndays_pvalues_rcp85.nc')
FWIpvalues_08 = xr.open_dataset(path_FWI_NorESM1_M_2+'FWI_ndays_pvalues_rcp85.nc')

# FWI trend (all models)
FWI = xr.concat((FWI_01['__xarray_dataarray_variable__'],FWI_02['__xarray_dataarray_variable__'],FWI_03['__xarray_dataarray_variable__'],FWI_04['__xarray_dataarray_variable__'],
                 FWI_05['fwi'],FWI_06['fwi'],FWI_07['fwi'],FWI_08['fwi']),dim='model')
# FWI trend pvalues (all models)
FWIpvalues = xr.concat((FWIpvalues_01['__xarray_dataarray_variable__'],FWIpvalues_02['__xarray_dataarray_variable__'],FWIpvalues_03['__xarray_dataarray_variable__'],FWIpvalues_04['__xarray_dataarray_variable__'],
                        FWIpvalues_05['fwi'],FWIpvalues_06['fwi'],FWIpvalues_07['fwi'],FWIpvalues_08['fwi']),dim='model')

# Model means
FWI_trend = FWI.mean(dim='model',skipna=True)
FWI_trend = xr.where(FWI_trend==0.,np.nan,FWI_trend)
FWI_min = FWI.min(dim='model',skipna=True)                             # Minimum trend between models
FWI_max = FWI.max(dim='model',skipna=True)                             # Maximum trend between models
FWIpvalues_max = FWIpvalues.max(dim='model',skipna=True)               # Maximum pvalues between models
FWIpvalues_max = xr.where(FWIpvalues_max<0.05,np.nan,FWIpvalues_max)   # Mask where maximum pvalues>=0.5
FWIpvalues_max = xr.where(FWIpvalues_max>=0.05,0.5,FWIpvalues_max)
FWItrend_sign = xr.full_like(FWIpvalues_max, fill_value=np.nan)        # Mask where trend sign is different between models
FWItrend_sign = xr.where(((FWI_min<0.) & (FWI_max>0.)),0.5,FWItrend_sign)
del FWI_min, FWI_max, FWI_01, FWI_02, FWI_03, FWI_04, FWI_05, FWI_06, FWI_07, FWI_08, FWIpvalues_01, FWIpvalues_02, FWIpvalues_03, FWIpvalues_04, FWIpvalues_05, FWIpvalues_06, FWIpvalues_07, FWIpvalues_08, FWI, FWIpvalues

########## Kindex ##########
path_Ki_EC_EARTH='../4-Projections_firestorm-risk/2-Kindex/COSMO-crCLIM-v1-1/EC-EARTH/'
path_Ki_HadGEM2_ES='../4-Projections_firestorm-risk/2-Kindex/COSMO-crCLIM-v1-1/HadGEM2-ES/'
path_Ki_MPI_ESM_LR='../4-Projections_firestorm-risk/2-Kindex/COSMO-crCLIM-v1-1/MPI-ESM-LR/'
path_Ki_NorESM1_M='../4-Projections_firestorm-risk/2-Kindex/COSMO-crCLIM-v1-1/NorESM1-M/'
path_Ki_EC_EARTH_2='../4-Projections_firestorm-risk/2-Kindex/RCA4/EC-EARTH/'
path_Ki_HadGEM2_ES_2='../4-Projections_firestorm-risk/2-Kindex/RCA4/HadGEM2-ES/'
path_Ki_MPI_ESM_LR_2='../4-Projections_firestorm-risk/2-Kindex/RCA4/MPI-ESM-LR/'
path_Ki_NorESM1_M_2='../4-Projections_firestorm-risk/2-Kindex/RCA4/NorESM1-M/'

# Read files
Ki_01 = xr.open_dataset(path_Ki_EC_EARTH+'Ki_ndays_trend_rcp85.nc')
Ki_02 = xr.open_dataset(path_Ki_HadGEM2_ES+'Ki_ndays_trend_rcp85.nc')
Ki_03 = xr.open_dataset(path_Ki_MPI_ESM_LR+'Ki_ndays_trend_rcp85.nc')
Ki_04 = xr.open_dataset(path_Ki_NorESM1_M+'Ki_ndays_trend_rcp85.nc')
Ki_05 = xr.open_dataset(path_Ki_EC_EARTH_2+'Ki_ndays_trend_rcp85.nc')
Ki_06 = xr.open_dataset(path_Ki_HadGEM2_ES_2+'Ki_ndays_trend_rcp85.nc')
Ki_07 = xr.open_dataset(path_Ki_MPI_ESM_LR_2+'Ki_ndays_trend_rcp85.nc')
Ki_08 = xr.open_dataset(path_Ki_NorESM1_M_2+'Ki_ndays_trend_rcp85.nc')
Kipvalues_01 = xr.open_dataset(path_Ki_EC_EARTH+'Ki_ndays_pvalues_rcp85.nc')
Kipvalues_02 = xr.open_dataset(path_Ki_HadGEM2_ES+'Ki_ndays_pvalues_rcp85.nc')
Kipvalues_03 = xr.open_dataset(path_Ki_MPI_ESM_LR+'Ki_ndays_pvalues_rcp85.nc')
Kipvalues_04 = xr.open_dataset(path_Ki_NorESM1_M+'Ki_ndays_pvalues_rcp85.nc')
Kipvalues_05 = xr.open_dataset(path_Ki_EC_EARTH_2+'Ki_ndays_pvalues_rcp85.nc')
Kipvalues_06 = xr.open_dataset(path_Ki_HadGEM2_ES_2+'Ki_ndays_pvalues_rcp85.nc')
Kipvalues_07 = xr.open_dataset(path_Ki_MPI_ESM_LR_2+'Ki_ndays_pvalues_rcp85.nc')
Kipvalues_08 = xr.open_dataset(path_Ki_NorESM1_M_2+'Ki_ndays_pvalues_rcp85.nc')

# Kindex trend (all models)
Ki = xr.concat((Ki_01['__xarray_dataarray_variable__'],Ki_02['__xarray_dataarray_variable__'],Ki_03['__xarray_dataarray_variable__'],Ki_04['__xarray_dataarray_variable__'],
                Ki_05['__xarray_dataarray_variable__'],Ki_06['__xarray_dataarray_variable__'],Ki_07['__xarray_dataarray_variable__'],Ki_08['__xarray_dataarray_variable__']),dim='model')
# Kindex trend pvalues (all models)
Kipvalues = xr.concat((Kipvalues_01['__xarray_dataarray_variable__'],Kipvalues_02['__xarray_dataarray_variable__'],Kipvalues_03['__xarray_dataarray_variable__'],Kipvalues_04['__xarray_dataarray_variable__'],
                       Kipvalues_05['__xarray_dataarray_variable__'],Kipvalues_06['__xarray_dataarray_variable__'],Kipvalues_07['__xarray_dataarray_variable__'],Kipvalues_08['__xarray_dataarray_variable__']),dim='model')

# Model means
Ki_trend = Ki.mean(dim='model',skipna=True)
Ki_trend = xr.where(Ki_trend==0.,np.nan,Ki_trend)
Ki_min = Ki.min(dim='model',skipna=True)                             # Minimum trend between models
Ki_max = Ki.max(dim='model',skipna=True)                             # Maximum trend between models
Kipvalues_max = Kipvalues.max(dim='model',skipna=True)               # Maximum pvalues between models
Kipvalues_max = xr.where(Kipvalues_max<0.05,np.nan,Kipvalues_max)    # Mask where maximum pvalues>=0.5
Kipvalues_max = xr.where(Kipvalues_max>=0.05,0.5,Kipvalues_max)
Kitrend_sign = xr.full_like(Kipvalues_max, fill_value=np.nan)        # Mask where trend sign is different between models
Kitrend_sign = xr.where(((Ki_min<0.) & (Ki_max>0.)),0.5,Kitrend_sign)
del Ki_min, Ki_max, Ki_01, Ki_02, Ki_03, Ki_04, Ki_05, Ki_06, Ki_07, Ki_08, Kipvalues_01, Kipvalues_02, Kipvalues_03, Kipvalues_04, Kipvalues_05, Kipvalues_06, Kipvalues_07, Kipvalues_08, Ki, Kipvalues

########## FWI and Kindex ##########
path_FWIKi_EC_EARTH='../4-Projections_firestorm-risk/3-Firestorm-risk/COSMO-crCLIM-v1-1/EC-EARTH/'
path_FWIKi_HadGEM2_ES='../4-Projections_firestorm-risk/3-Firestorm-risk/COSMO-crCLIM-v1-1/HadGEM2-ES/'
path_FWIKi_MPI_ESM_LR='../4-Projections_firestorm-risk/3-Firestorm-risk/COSMO-crCLIM-v1-1/MPI-ESM-LR/'
path_FWIKi_NorESM1_M='../4-Projections_firestorm-risk/3-Firestorm-risk/COSMO-crCLIM-v1-1/NorESM1-M/'
path_FWIKi_EC_EARTH_2='../4-Projections_firestorm-risk/3-Firestorm-risk/RCA4/EC-EARTH/'
path_FWIKi_HadGEM2_ES_2='../4-Projections_firestorm-risk/3-Firestorm-risk/RCA4/HadGEM2-ES/'
path_FWIKi_MPI_ESM_LR_2='../4-Projections_firestorm-risk/3-Firestorm-risk/RCA4/MPI-ESM-LR/'
path_FWIKi_NorESM1_M_2='../4-Projections_firestorm-risk/3-Firestorm-risk/RCA4/NorESM1-M/'

# Read files
FWIKi_01 = xr.open_dataset(path_FWIKi_EC_EARTH+'FWI-Ki_ndays_trend_rcp85.nc')
FWIKi_02 = xr.open_dataset(path_FWIKi_HadGEM2_ES+'FWI-Ki_ndays_trend_rcp85.nc')
FWIKi_03 = xr.open_dataset(path_FWIKi_MPI_ESM_LR+'FWI-Ki_ndays_trend_rcp85.nc')
FWIKi_04 = xr.open_dataset(path_FWIKi_NorESM1_M+'FWI-Ki_ndays_trend_rcp85.nc')
FWIKi_05 = xr.open_dataset(path_FWIKi_EC_EARTH_2+'FWI-Ki_ndays_trend_rcp85.nc')
FWIKi_06 = xr.open_dataset(path_FWIKi_HadGEM2_ES_2+'FWI-Ki_ndays_trend_rcp85.nc')
FWIKi_07 = xr.open_dataset(path_FWIKi_MPI_ESM_LR_2+'FWI-Ki_ndays_trend_rcp85.nc')
FWIKi_08 = xr.open_dataset(path_FWIKi_NorESM1_M_2+'FWI-Ki_ndays_trend_rcp85.nc')
FWIKipvalues_01 = xr.open_dataset(path_FWIKi_EC_EARTH+'FWI-Ki_ndays_pvalues_rcp85.nc')
FWIKipvalues_02 = xr.open_dataset(path_FWIKi_HadGEM2_ES+'FWI-Ki_ndays_pvalues_rcp85.nc')
FWIKipvalues_03 = xr.open_dataset(path_FWIKi_MPI_ESM_LR+'FWI-Ki_ndays_pvalues_rcp85.nc')
FWIKipvalues_04 = xr.open_dataset(path_FWIKi_NorESM1_M+'FWI-Ki_ndays_pvalues_rcp85.nc')
FWIKipvalues_05 = xr.open_dataset(path_FWIKi_EC_EARTH_2+'FWI-Ki_ndays_pvalues_rcp85.nc')
FWIKipvalues_06 = xr.open_dataset(path_FWIKi_HadGEM2_ES_2+'FWI-Ki_ndays_pvalues_rcp85.nc')
FWIKipvalues_07 = xr.open_dataset(path_FWIKi_MPI_ESM_LR_2+'FWI-Ki_ndays_pvalues_rcp85.nc')
FWIKipvalues_08 = xr.open_dataset(path_FWIKi_NorESM1_M_2+'FWI-Ki_ndays_pvalues_rcp85.nc')

# FWI and Kindex trend (all models)
FWIKi = xr.concat((FWIKi_01['__xarray_dataarray_variable__'],FWIKi_02['__xarray_dataarray_variable__'],FWIKi_03['__xarray_dataarray_variable__'],FWIKi_04['__xarray_dataarray_variable__'],
                   FWIKi_05['__xarray_dataarray_variable__'],FWIKi_06['__xarray_dataarray_variable__'],FWIKi_07['__xarray_dataarray_variable__'],FWIKi_08['__xarray_dataarray_variable__']),dim='model')
# FWI and Kindex trend pvalues (all models)
FWIKipvalues = xr.concat((FWIKipvalues_01['__xarray_dataarray_variable__'],FWIKipvalues_02['__xarray_dataarray_variable__'],FWIKipvalues_03['__xarray_dataarray_variable__'],FWIKipvalues_04['__xarray_dataarray_variable__'],
                          FWIKipvalues_05['__xarray_dataarray_variable__'],FWIKipvalues_06['__xarray_dataarray_variable__'],FWIKipvalues_07['__xarray_dataarray_variable__'],FWIKipvalues_08['__xarray_dataarray_variable__']),dim='model')

# Model means
FWIKi_trend = FWIKi.mean(dim='model',skipna=True)
FWIKi_trend = xr.where(FWIKi_trend==0.,np.nan,FWIKi_trend)
FWIKi_min = FWIKi.min(dim='model',skipna=True)                             # Minimum trend between models
FWIKi_max = FWIKi.max(dim='model',skipna=True)                             # Maximum trend between models
FWIKipvalues_max = FWIKipvalues.max(dim='model',skipna=True)               # Maximum pvalues between models
FWIKipvalues_max = xr.where(FWIKipvalues_max<0.05,np.nan,FWIKipvalues_max) # Mask where maximum pvalues>=0.5
FWIKipvalues_max = xr.where(FWIKipvalues_max>=0.05,0.5,FWIKipvalues_max)
FWIKitrend_sign = xr.full_like(FWIKipvalues_max, fill_value=np.nan)        # Mask where trend sign is different between models
FWIKitrend_sign = xr.where(((FWIKi_min<0.) & (FWIKi_max>0.)),0.5,FWIKitrend_sign)
del FWIKi_min, FWIKi_max, FWIKi_01, FWIKi_02, FWIKi_03, FWIKi_04, FWIKi_05, FWIKi_06, FWIKi_07, FWIKi_08, FWIKipvalues_01, FWIKipvalues_02, FWIKipvalues_03, FWIKipvalues_04, FWIKipvalues_05, FWIKipvalues_06, FWIKipvalues_07, FWIKipvalues_08, FWIKi, FWIKipvalues

########## Country time-series ##########
path_countries='./Utils/'
def ndayscalc(count):
      COUNT = xr.open_dataset(path_countries+count+'.nc')['Band1']
      mask = COUNT.interp(lat=FWI_trend.lat, lon=FWI_trend.lon, method='linear') # Interpolation

      ########## FWI and Kindex ##########
      Nd_01 = xr.open_dataset(path_FWIKi_EC_EARTH+'FWI-Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_02 = xr.open_dataset(path_FWIKi_HadGEM2_ES+'FWI-Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_03 = xr.open_dataset(path_FWIKi_MPI_ESM_LR+'FWI-Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_04 = xr.open_dataset(path_FWIKi_NorESM1_M+'FWI-Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_05 = xr.open_dataset(path_FWIKi_EC_EARTH_2+'FWI-Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_06 = xr.open_dataset(path_FWIKi_HadGEM2_ES_2+'FWI-Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_07 = xr.open_dataset(path_FWIKi_MPI_ESM_LR_2+'FWI-Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_08 = xr.open_dataset(path_FWIKi_NorESM1_M_2+'FWI-Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))

      # Model mean and standard deviation
      Nd = xr.concat((Nd_01,Nd_02,Nd_03,Nd_04,Nd_05,Nd_06,Nd_07,Nd_08),dim='model')
      Nd_mean = Nd.mean(dim='model',skipna=True)
      Nd_std = Nd.std(dim='model',skipna=True)
      del Nd, Nd_01, Nd_02, Nd_03, Nd_04, Nd_05, Nd_06, Nd_07, Nd_08

      ########## FWI ##########
      Nd_01 = xr.open_dataset(path_FWI_EC_EARTH+'FWI_ndays_rcp85.nc')['fwi'].where(mask==1).mean(axis=(1,2))
      Nd_02 = xr.open_dataset(path_FWI_HadGEM2_ES+'FWI_ndays_rcp85.nc')['fwi'].where(mask==1).mean(axis=(1,2))
      Nd_03 = xr.open_dataset(path_FWI_MPI_ESM_LR+'FWI_ndays_rcp85.nc')['fwi'].where(mask==1).mean(axis=(1,2))
      Nd_04 = xr.open_dataset(path_FWI_NorESM1_M+'FWI_ndays_rcp85.nc')['fwi'].where(mask==1).mean(axis=(1,2))
      Nd_05 = xr.open_dataset(path_FWI_EC_EARTH_2+'FWI_ndays_rcp85.nc')['fwi'].where(mask==1).mean(axis=(1,2))
      Nd_06 = xr.open_dataset(path_FWI_HadGEM2_ES_2+'FWI_ndays_rcp85.nc')['fwi'].where(mask==1).mean(axis=(1,2))
      Nd_07 = xr.open_dataset(path_FWI_MPI_ESM_LR_2+'FWI_ndays_rcp85.nc')['fwi'].where(mask==1).mean(axis=(1,2))
      Nd_08 = xr.open_dataset(path_FWI_NorESM1_M_2+'FWI_ndays_rcp85.nc')['fwi'].where(mask==1).mean(axis=(1,2))

      # Model mean and standard deviation
      Nd = xr.concat((Nd_01,Nd_02,Nd_03,Nd_04,Nd_05,Nd_06,Nd_07,Nd_08),dim='model')
      Nd_fwi_mean = Nd.mean(dim='model',skipna=True)
      Nd_fwi_std = Nd.std(dim='model',skipna=True)
      del Nd, Nd_01, Nd_02, Nd_03, Nd_04, Nd_05, Nd_06, Nd_07, Nd_08

      ########## Kindex ##########
      Nd_01 = xr.open_dataset(path_Ki_EC_EARTH+'Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_02 = xr.open_dataset(path_Ki_HadGEM2_ES+'Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_03 = xr.open_dataset(path_Ki_MPI_ESM_LR+'Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_04 = xr.open_dataset(path_Ki_NorESM1_M+'Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_05 = xr.open_dataset(path_Ki_EC_EARTH_2+'Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_06 = xr.open_dataset(path_Ki_HadGEM2_ES_2+'Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_07 = xr.open_dataset(path_Ki_MPI_ESM_LR_2+'Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))
      Nd_08 = xr.open_dataset(path_Ki_NorESM1_M_2+'Ki_ndays_rcp85.nc')['__xarray_dataarray_variable__'].where(mask==1).mean(axis=(1,2))

      # Model mean and standard deviation
      Nd = xr.concat((Nd_01,Nd_02,Nd_03,Nd_04,Nd_05,Nd_06,Nd_07,Nd_08),dim='model')
      Nd_ki_mean = Nd.mean(dim='model',skipna=True)
      Nd_ki_std = Nd.std(dim='model',skipna=True)
      del Nd, Nd_01, Nd_02, Nd_03, Nd_04, Nd_05, Nd_06, Nd_07, Nd_08

      return Nd_mean, Nd_std, Nd_fwi_mean, Nd_fwi_std, Nd_ki_mean, Nd_ki_std

Nd_mean_PT, Nd_std_PT, Nd_fwi_mean_PT, Nd_fwi_std_PT, Nd_ki_mean_PT, Nd_ki_std_PT = ndayscalc('Portugal')
Nd_mean_SP, Nd_std_SP, Nd_fwi_mean_SP, Nd_fwi_std_SP, Nd_ki_mean_SP, Nd_ki_std_SP = ndayscalc('Spain')
Nd_mean_FR, Nd_std_FR, Nd_fwi_mean_FR, Nd_fwi_std_FR, Nd_ki_mean_FR, Nd_ki_std_FR = ndayscalc('France')
Nd_mean_IT, Nd_std_IT, Nd_fwi_mean_IT, Nd_fwi_std_IT, Nd_ki_mean_IT, Nd_ki_std_IT = ndayscalc('Italy')
Nd_mean_GR, Nd_std_GR, Nd_fwi_mean_GR, Nd_fwi_std_GR, Nd_ki_mean_GR, Nd_ki_std_GR = ndayscalc('Greece')
Nd_mean_TU, Nd_std_TU, Nd_fwi_mean_TU, Nd_fwi_std_TU, Nd_ki_mean_TU, Nd_ki_std_TU = ndayscalc('Turkey')
dates = Nd_mean_PT.time['time.year'].time

# Map definition
m = Basemap(projection='tmerc', resolution='h',
                        llcrnrlat=30., urcrnrlat=65.,
                        llcrnrlon=-10., urcrnrlon=62.,
                        lat_0=47., lon_0=15.)
lats = FWI_trend['lat']
lons = FWI_trend['lon']
x, y = m(lons, lats)

# Plot of the FWI trend
fig = plt.figure(figsize=(5,7))
gs = fig.add_gridspec(6,2)
ax1 = fig.add_subplot(gs[0:2,0])
m.drawcoastlines(linewidth=0.4,zorder=3)
m.drawcountries(linewidth=0.3,zorder=3)
m.drawparallels(np.arange(0,80,5),linewidth=0.3,labels=[1,0,0,0],zorder=3)
m.drawmeridians(np.arange(-60,80,10),linewidth=0.3,labels=[0,0,0,1],zorder=3)
m.drawlsmask(land_color='none',ocean_color='white',lakes=True,zorder=2)
clevs =  np.arange(-5.,5.5,0.5)
cs2 = m.contourf(x, y, FWI_trend*10., clevs, cmap=cmaps.ViBlGrWhYeOrRe, alpha=0.95, extend='both',zorder=1)
cs3 = m.contourf(x, y, FWIpvalues_max, np.arange(0.,1.1,0.1), cmap=get_cmap("binary"), alpha=1., zorder=1)
cs4 = m.contourf(x, y, FWItrend_sign, np.arange(0.,1.1,0.1), cmap=get_cmap('binary'), alpha=1., zorder=1)
plt.title('FWI trend ($days/decade$)',loc='left',weight='bold',size=5,pad=4.)
ax1.text(-0.1, 1.05, 'a',family='sans-serif',weight='bold',size=7, horizontalalignment='left', verticalalignment='top',transform=ax1.transAxes)
axins = inset_axes(ax1,
                   width="3%",  
                   height="100%",  
                   loc='lower left',
                   bbox_to_anchor=(1.01,0., 1, 1),
                   bbox_transform=ax1.transAxes,
                   borderpad=0,
                   )
cbar = fig.colorbar(cs2,orientation="vertical",cax=axins,ticks=np.arange(-5.,6.,1.))
plt.gcf()

# Plot of the Kindex trend
ax2 = fig.add_subplot(gs[2:4,0])
m.drawcoastlines(linewidth=0.4,zorder=3)
m.drawcountries(linewidth=0.3,zorder=3)
m.drawparallels(np.arange(0,80,5),linewidth=0.3,labels=[1,0,0,0],zorder=3)
m.drawmeridians(np.arange(-60,80,10),linewidth=0.3,labels=[0,0,0,1],zorder=3)
m.drawlsmask(land_color='none',ocean_color='white',lakes=True,zorder=2)
clevs =  np.arange(-5.,5.5,0.5)
cs2 = m.contourf(x, y, Ki_trend*10., clevs, cmap=cmaps.ViBlGrWhYeOrRe, alpha=0.95, extend='both',zorder=1)
cs3 = m.contourf(x, y, Kipvalues_max, np.arange(0.,1.1,0.1), cmap=get_cmap("binary"), alpha=1., zorder=1)
cs4 = m.contourf(x, y, Kitrend_sign, np.arange(0.,1.1,0.1), cmap=get_cmap("binary"), alpha=1., zorder=1)
plt.title('K-index trend ($days/decade$)',loc='left',weight='bold',size=5,pad=4.)
ax2.text(-0.1, 1.05, 'b',family='sans-serif',weight='bold',size=7, horizontalalignment='left', verticalalignment='top',transform=ax2.transAxes)
axins = inset_axes(ax2,
                   width="3%",  
                   height="100%",  
                   loc='lower left',
                   bbox_to_anchor=(1.01,0., 1, 1),
                   bbox_transform=ax2.transAxes,
                   borderpad=0,
                   )
cbar = fig.colorbar(cs2,orientation="vertical",cax=axins,ticks= np.arange(-5.,6.,1.))

# Plot of the firestorm risk trend
ax3 = fig.add_subplot(gs[4:6,0])
m.drawcoastlines(linewidth=0.4,zorder=3)
m.drawcountries(linewidth=0.3,zorder=3)
m.drawparallels(np.arange(0,80,5),linewidth=0.3,labels=[1,0,0,0],zorder=3)
m.drawmeridians(np.arange(-60,80,10),linewidth=0.3,labels=[0,0,0,1],zorder=3)
m.drawlsmask(land_color='none',ocean_color='white',lakes=True,zorder=2)
clevs =  np.arange(-5.,5.5,0.5)
cs2 = m.contourf(x, y, FWIKi_trend*10., clevs, cmap=cmaps.ViBlGrWhYeOrRe, alpha=0.95, extend='both',zorder=1)
cs3 = m.contourf(x, y, FWIKipvalues_max, np.arange(0.,1.1,0.1), cmap=get_cmap("binary"), alpha=1., zorder=1)
cs4 = m.contourf(x, y, FWIKitrend_sign, np.arange(0.,1.1,0.1), cmap=get_cmap("binary"), alpha=1., zorder=1)
plt.title('Atmospheric firestorm risk trend ($days/decade$)',loc='left',weight='bold',size=5,pad=4.)
ax3.text(-0.1, 1.05, 'c',family='sans-serif',weight='bold',size=7, horizontalalignment='left', verticalalignment='top',transform=ax3.transAxes)
axins = inset_axes(ax3,
                   width="3%",  
                   height="100%",  
                   loc='lower left',
                   bbox_to_anchor=(1.01,0., 1, 1),
                   bbox_transform=ax3.transAxes,
                   borderpad=0,
                   )
cbar = fig.colorbar(cs2,orientation="vertical",cax=axins,ticks= np.arange(-5.,6.,1.))
plt.gcf()

# Plot of the Portugal time-series
ax4 = fig.add_subplot(gs[0,1])
ax4.fill_between(dates, Nd_mean_PT-Nd_std_PT, Nd_mean_PT+Nd_std_PT, facecolor='k', alpha=0.5, interpolate=True)
ax4.plot(dates,Nd_mean_PT,'k',linewidth=0.7)
ax4.fill_between(dates, Nd_fwi_mean_PT-Nd_fwi_std_PT, Nd_fwi_mean_PT+Nd_fwi_std_PT, facecolor='r', alpha=0.5, interpolate=True)
ax4.plot(dates,Nd_fwi_mean_PT,'r',linewidth=0.7)
ax4.fill_between(dates, Nd_ki_mean_PT-Nd_ki_std_PT, Nd_ki_mean_PT+Nd_ki_std_PT, facecolor='b', alpha=0.5, interpolate=True)
ax4.plot(dates,Nd_ki_mean_PT,'b',linewidth=0.7)
ax4.set_ylim([0,150])
ax4.xaxis.set_minor_locator(matplotlib.dates.YearLocator(10))
ax4.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))
ax4.xaxis.set_major_locator(matplotlib.dates.YearLocator(20))
plt.title('Portugal',loc='left',weight='bold',size=5,pad=1.)
ax4.margins(x=0, y=0)
ax4.tick_params(axis='both', which='both', length=4., pad=1., labelsize=5, width=0.5)
ax4.grid(True,which='both',linestyle='-',linewidth=0.5)
plt.gcf()

# Plot of the Spain time-series
ax5 = fig.add_subplot(gs[1,1])
ax5.fill_between(dates, Nd_mean_SP-Nd_std_SP, Nd_mean_SP+Nd_std_SP, facecolor='k', alpha=0.5, interpolate=True)
ax5.plot(dates,Nd_mean_SP,'k',linewidth=0.7)
ax5.fill_between(dates, Nd_fwi_mean_SP-Nd_fwi_std_SP, Nd_fwi_mean_SP+Nd_fwi_std_SP, facecolor='r', alpha=0.5, interpolate=True)
ax5.plot(dates,Nd_fwi_mean_SP,'r',linewidth=0.7)
ax5.fill_between(dates, Nd_ki_mean_SP-Nd_ki_std_SP, Nd_ki_mean_SP+Nd_ki_std_SP, facecolor='b', alpha=0.5, interpolate=True)
ax5.plot(dates,Nd_ki_mean_SP,'b',linewidth=0.7)
ax5.set_ylim([0,150])
ax5.xaxis.set_minor_locator(matplotlib.dates.YearLocator(10))
ax5.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))
ax5.xaxis.set_major_locator(matplotlib.dates.YearLocator(20))
plt.title('Spain',loc='left',weight='bold',size=5,pad=1.)
ax5.margins(x=0, y=0)
ax5.tick_params(axis='both', which='both', length=4., pad=1., labelsize=5, width=0.5)
ax5.grid(True,which='both',linestyle='-',linewidth=0.5)
plt.gcf()

# Plot of the France time-series
ax6 = fig.add_subplot(gs[2,1])
ax6.fill_between(dates, Nd_mean_FR-Nd_std_FR, Nd_mean_FR+Nd_std_FR, facecolor='k', alpha=0.5, interpolate=True)
ax6.plot(dates,Nd_mean_FR,'k',linewidth=0.7)
ax6.fill_between(dates, Nd_fwi_mean_FR-Nd_fwi_std_FR, Nd_fwi_mean_FR+Nd_fwi_std_FR, facecolor='r', alpha=0.5, interpolate=True)
ax6.plot(dates,Nd_fwi_mean_FR,'r',linewidth=0.7)
ax6.fill_between(dates, Nd_ki_mean_FR-Nd_ki_std_FR, Nd_ki_mean_FR+Nd_ki_std_FR, facecolor='b', alpha=0.5, interpolate=True)
ax6.plot(dates,Nd_ki_mean_FR,'b',linewidth=0.7)
ax6.set_ylim([0,150])
ax6.xaxis.set_minor_locator(matplotlib.dates.YearLocator(10))
ax6.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))
ax6.xaxis.set_major_locator(matplotlib.dates.YearLocator(20))
plt.title('France',loc='left',weight='bold',size=5,pad=1.)
ax6.margins(x=0, y=0)
ax6.tick_params(axis='both', which='both', length=4., pad=1., labelsize=5, width=0.5)
ax6.grid(True,which='both',linestyle='-',linewidth=0.5)
plt.gcf()

# Plot of the Italy time-series
ax7 = fig.add_subplot(gs[3,1])
ax7.fill_between(dates, Nd_mean_IT-Nd_std_IT, Nd_mean_IT+Nd_std_IT, facecolor='k', alpha=0.5, interpolate=True)
ax7.plot(dates,Nd_mean_IT,'k',linewidth=0.7)
ax7.fill_between(dates, Nd_fwi_mean_IT-Nd_fwi_std_IT, Nd_fwi_mean_IT+Nd_fwi_std_IT, facecolor='r', alpha=0.5, interpolate=True)
ax7.plot(dates,Nd_fwi_mean_IT,'r',linewidth=0.7)
ax7.fill_between(dates, Nd_ki_mean_IT-Nd_ki_std_IT, Nd_ki_mean_IT+Nd_ki_std_IT, facecolor='b', alpha=0.5, interpolate=True)
ax7.plot(dates,Nd_ki_mean_IT,'b',linewidth=0.7)
ax7.set_ylim([0,150])
ax7.xaxis.set_minor_locator(matplotlib.dates.YearLocator(10))
ax7.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))
ax7.xaxis.set_major_locator(matplotlib.dates.YearLocator(20))
plt.title('Italy',loc='left',weight='bold',size=5,pad=1.)
ax7.margins(x=0, y=0)
ax7.tick_params(axis='both', which='both', length=4., pad=1., labelsize=5, width=0.5)
ax7.grid(True,which='both',linestyle='-',linewidth=0.5)
plt.gcf()

# Plot of the Greece time-series
ax8 = fig.add_subplot(gs[4,1])
ax8.fill_between(dates, Nd_mean_GR-Nd_std_GR, Nd_mean_GR+Nd_std_GR, facecolor='k', alpha=0.5, interpolate=True)
ax8.plot(dates,Nd_mean_GR,'k',linewidth=0.7)
ax8.fill_between(dates, Nd_fwi_mean_GR-Nd_fwi_std_GR, Nd_fwi_mean_GR+Nd_fwi_std_GR, facecolor='r', alpha=0.5, interpolate=True)
ax8.plot(dates,Nd_fwi_mean_GR,'r',linewidth=0.7)
ax8.fill_between(dates, Nd_ki_mean_GR-Nd_ki_std_GR, Nd_ki_mean_GR+Nd_ki_std_GR, facecolor='b', alpha=0.5, interpolate=True)
ax8.plot(dates,Nd_ki_mean_GR,'b',linewidth=0.7)
ax8.set_ylim([0,150])
ax8.xaxis.set_minor_locator(matplotlib.dates.YearLocator(10))
ax8.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))
ax8.xaxis.set_major_locator(matplotlib.dates.YearLocator(20))
plt.title('Greece',loc='left',weight='bold',size=5,pad=1.)
ax8.margins(x=0, y=0)
ax8.tick_params(axis='both', which='both', length=4., pad=1., labelsize=5, width=0.5)
ax8.grid(True,which='both',linestyle='-',linewidth=0.5)
plt.gcf()

# Plot of the Turkey time-series
ax9 = fig.add_subplot(gs[5,1])
ax9.fill_between(dates, Nd_mean_TU-Nd_std_TU, Nd_mean_TU+Nd_std_TU, facecolor='k', alpha=0.5, interpolate=True)
ax9.plot(dates,Nd_mean_TU,'k',linewidth=0.7)
ax9.fill_between(dates, Nd_fwi_mean_TU-Nd_fwi_std_TU, Nd_fwi_mean_TU+Nd_fwi_std_TU, facecolor='r', alpha=0.5, interpolate=True)
ax9.plot(dates,Nd_fwi_mean_TU,'r',linewidth=0.7)
ax9.fill_between(dates, Nd_ki_mean_TU-Nd_ki_std_TU, Nd_ki_mean_TU+Nd_ki_std_TU, facecolor='b', alpha=0.5, interpolate=True)
ax9.plot(dates,Nd_ki_mean_TU,'b',linewidth=0.7)
ax9.set_ylim([0,150])
ax9.xaxis.set_minor_locator(matplotlib.dates.YearLocator(10))
ax9.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%Y'))
ax9.xaxis.set_major_locator(matplotlib.dates.YearLocator(20))
plt.title('Turkey',loc='left',weight='bold',size=5,pad=1.)
ax9.margins(x=0, y=0)
ax9.tick_params(axis='both', which='both', length=3., pad=1., labelsize=5, width=0.5)
ax9.grid(True,which='both',linestyle='-',linewidth=0.5)
plt.gcf()

# Save figure
plt.subplots_adjust(top = 0.9, bottom=0.1, hspace=0.3, wspace=0.3)
fig.savefig('Figure4.png',dpi=300,bbox_inches='tight')   
fig.savefig('Figure4.pdf',bbox_inches='tight')   
